Let us set some global options for all code chunks in this document.

knitr::opts_chunk$set(
  message = FALSE,    # Disable messages printed by R code chunks
  warning = FALSE,    # Disable warnings printed by R code chunks
  echo = TRUE,        # Show R code within code chunks in output
  include = TRUE,     # Include both R code and its results in output
  eval = TRUE,       # Evaluate R code chunks
  cache = FALSE,       # Enable caching of R code chunks for faster rendering
  fig.align = "center",
  out.width = "100%",
  retina = 2,
  error = TRUE,
  collapse = TRUE
)
rm(list = ls())
set.seed(1982)

1 Preprocessing

Let us now load some required libraries.

# Load required libraries

# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")

library(INLA)
inla.setOption(num.threads = 7)
library(inlabru)
library(rSPDE)
library(MetricGraph)

library(plotly)
library(dplyr)

library(sf)

library(here)

Function standarize() below is later used to standardize the covariate SpeedLimit.

standardize <- function(x) {return((x - mean(x)) / sd(x))}

To keep track of the changes, we provide summaries of every new created object. Those summaries can be accessed by pressing the Show buttons below


We load the graph object sf_graph (which only contains weights) and the data (already graph-processed).

load(here("Graph_objects/graph_construction_25_04_2024partialtomtomwhichlonglatsf.RData"))
load(here("Data_files/data_day7142128_hour13_with_no_consecutive_zeros_partialtomtom_graph_25_04_2024_processed.RData"))
data_on_graph = data_on_graph %>% 
  dplyr::select(-datetime)
sf_graph$get_edge_lengths() %>% head() %>% capture.output() %>% grep("^Units:", ., value = TRUE)
## [1] "Units: [km]"
summary(sf_graph)
## A metric graph object with:
## 
## Vertices:
##   Total: 4017 
##   Degree 2: 1821;  Degree 3: 180;  Degree 4: 1402;  Degree 5: 94;  Degree 6: 367; 
##   Degree 7: 36;  Degree 8: 116;  Degree 12: 1; 
##   With incompatible directions:  0 
## 
## Edges: 
##   Total: 6827 
##   Lengths: 
##       Min: 0.002834658  ; Max: 0.2743201  ; Total: 311.1441 
##   Weights: 
##       Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour 
##   That are circles:  0 
## 
## Graph units: 
##   Vertices unit:  degrees  ; Lengths unit:  km 
## 
## Longitude and Latitude coordinates:  TRUE
##   Which spatial package:  sf 
##   CRS:  EPSG:4326
## 
## Some characteristics of the graph:
##   Connected: TRUE
##   Has loops: FALSE
##   Has multiple edges: TRUE
##   Is a tree: FALSE
##   Distance consistent: TRUE
##   Has Euclidean edges: FALSE
## 
## Computed quantities inside the graph: 
##   Laplacian:  FALSE  ; Geodesic distances:  TRUE 
##   Resistance distances:  FALSE  ; Finite element matrices:  FALSE 
## 
## Mesh: The graph has no mesh! 
## 
## Data: The graph has no data!
## 
## Tolerances: 
##   vertex-vertex:  0.001 
##   vertex-edge:  0.001 
##   edge-edge:  0
summary(data_on_graph)
##        ID           speed             day        .distance_to_graph
##  Min.   :5701   Min.   : 0.000   Min.   :1.000   Min.   :0.000000  
##  1st Qu.:6571   1st Qu.: 1.609   1st Qu.:2.000   1st Qu.:0.001655  
##  Median :6687   Median :14.484   Median :3.000   Median :0.003552  
##  Mean   :7106   Mean   :15.142   Mean   :2.556   Mean   :0.004469  
##  3rd Qu.:7281   3rd Qu.:24.140   3rd Qu.:4.000   3rd Qu.:0.006152  
##  Max.   :8969   Max.   :99.779   Max.   :4.000   Max.   :0.019991  
##   .edge_number  .distance_on_edge    .group             .coord_x     
##  Min.   :   1   Min.   :0.0000    Length:39574       Min.   :-122.4  
##  1st Qu.:1413   1st Qu.:0.2726    Class :character   1st Qu.:-122.4  
##  Median :2879   Median :0.5149    Mode  :character   Median :-122.4  
##  Mean   :3090   Mean   :0.5089                       Mean   :-122.4  
##  3rd Qu.:4753   3rd Qu.:0.7501                       3rd Qu.:-122.4  
##  Max.   :6826   Max.   :1.0000                       Max.   :-122.4  
##     .coord_y    
##  Min.   :37.77  
##  1st Qu.:37.78  
##  Median :37.79  
##  Mean   :37.79  
##  3rd Qu.:37.79  
##  Max.   :37.81

The following commands remove zero speed observations that are 1m away from the graph, and after that, they remove any speed observations that are 3m away from the graph.

to_remove = data_on_graph %>%
  filter(speed == 0, .distance_to_graph > 0.001) 

data_on_graph = setdiff(data_on_graph, to_remove) %>% 
  filter(.distance_to_graph <= 0.003)
summary(to_remove)
##        ID           speed        day        .distance_to_graph  .edge_number 
##  Min.   :5701   Min.   :0   Min.   :1.000   Min.   :0.001001   Min.   :   1  
##  1st Qu.:6574   1st Qu.:0   1st Qu.:2.000   1st Qu.:0.002757   1st Qu.:1376  
##  Median :6688   Median :0   Median :3.000   Median :0.004563   Median :2788  
##  Mean   :7078   Mean   :0   Mean   :2.556   Mean   :0.005453   Mean   :3031  
##  3rd Qu.:7277   3rd Qu.:0   3rd Qu.:4.000   3rd Qu.:0.007129   3rd Qu.:4715  
##  Max.   :8969   Max.   :0   Max.   :4.000   Max.   :0.019988   Max.   :6812  
##  .distance_on_edge    .group             .coord_x         .coord_y    
##  Min.   :0.0000    Length:8597        Min.   :-122.4   Min.   :37.77  
##  1st Qu.:0.3078    Class :character   1st Qu.:-122.4   1st Qu.:37.78  
##  Median :0.5301    Mode  :character   Median :-122.4   Median :37.79  
##  Mean   :0.5177                       Mean   :-122.4   Mean   :37.79  
##  3rd Qu.:0.7362                       3rd Qu.:-122.4   3rd Qu.:37.79  
##  Max.   :1.0000                       Max.   :-122.4   Max.   :37.81
summary(data_on_graph)
##        ID           speed             day       .distance_to_graph 
##  Min.   :5701   Min.   : 0.000   Min.   :1.00   Min.   :0.0000000  
##  1st Qu.:6565   1st Qu.: 9.656   1st Qu.:2.00   1st Qu.:0.0005835  
##  Median :6683   Median :19.312   Median :3.00   Median :0.0012618  
##  Mean   :7115   Mean   :19.385   Mean   :2.55   Mean   :0.0013511  
##  3rd Qu.:7286   3rd Qu.:27.359   3rd Qu.:4.00   3rd Qu.:0.0021077  
##  Max.   :8969   Max.   :99.779   Max.   :4.00   Max.   :0.0029999  
##   .edge_number  .distance_on_edge    .group             .coord_x     
##  Min.   :   1   Min.   :0.0000    Length:14534       Min.   :-122.4  
##  1st Qu.:1376   1st Qu.:0.2552    Class :character   1st Qu.:-122.4  
##  Median :2956   Median :0.5139    Mode  :character   Median :-122.4  
##  Mean   :3130   Mean   :0.5099                       Mean   :-122.4  
##  3rd Qu.:4803   3rd Qu.:0.7697                       3rd Qu.:-122.4  
##  Max.   :6817   Max.   :1.0000                       Max.   :-122.4  
##     .coord_y    
##  Min.   :37.77  
##  1st Qu.:37.78  
##  Median :37.79  
##  Mean   :37.79  
##  3rd Qu.:37.79  
##  Max.   :37.81

We add data to the graph.

sf_graph$add_observations(data = data_on_graph, 
                          group = "day", 
                          normalized = TRUE, 
                          clear_obs = TRUE)
sf_graph$get_data()
## # A tibble: 14,534 × 9
##       ID speed   day .distance_to_graph .coord_x .coord_y .edge_number
##    <int> <dbl> <dbl>              <dbl>    <dbl>    <dbl>        <dbl>
##  1  6666   0       1           0.00100     -122.     37.8            1
##  2  8941  12.9     1           0.00230     -122.     37.8            4
##  3  8768  24.1     1           0.00233     -122.     37.8            6
##  4  8929  32.2     1           0.00151     -122.     37.8            6
##  5  8965   0       1           0.000647    -122.     37.8            9
##  6  8965  19.3     1           0.00170     -122.     37.8           14
##  7  8954  22.5     1           0.00103     -122.     37.8           14
##  8  8774  19.3     1           0.00279     -122.     37.8           14
##  9  6655  30.6     1           0.00208     -122.     37.8           18
## 10  6677  14.5     1           0.000436    -122.     37.8           20
## # ℹ 14,524 more rows
## # ℹ 2 more variables: .distance_on_edge <dbl>, .group <chr>
summary(sf_graph)
## A metric graph object with:
## 
## Vertices:
##   Total: 4017 
##   Degree 2: 1821;  Degree 3: 180;  Degree 4: 1402;  Degree 5: 94;  Degree 6: 367; 
##   Degree 7: 36;  Degree 8: 116;  Degree 12: 1; 
##   With incompatible directions:  0 
## 
## Edges: 
##   Total: 6827 
##   Lengths: 
##       Min: 0.002834658  ; Max: 0.2743201  ; Total: 311.1441 
##   Weights: 
##       Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour 
##   That are circles:  0 
## 
## Graph units: 
##   Vertices unit:  degrees  ; Lengths unit:  km 
## 
## Longitude and Latitude coordinates:  TRUE
##   Which spatial package:  sf 
##   CRS:  EPSG:4326
## 
## Some characteristics of the graph:
##   Connected: TRUE
##   Has loops: FALSE
##   Has multiple edges: TRUE
##   Is a tree: FALSE
##   Distance consistent: TRUE
##   Has Euclidean edges: FALSE
## 
## Computed quantities inside the graph: 
##   Laplacian:  FALSE  ; Geodesic distances:  TRUE 
##   Resistance distances:  FALSE  ; Finite element matrices:  FALSE 
## 
## Mesh: The graph has no mesh! 
## 
## Data: 
##   Columns:  ID speed day 
##   Groups:  .group 
## 
## Tolerances: 
##   vertex-vertex:  0.001 
##   vertex-edge:  0.001 
##   edge-edge:  0

We get the values of the weights at data locations. This essentially gives us covariates from the weights.

sf_graph$edgeweight_to_data(data_loc = TRUE)
sf_graph$get_data()
## # A tibble: 57,108 × 54
##       ID speed   day .distance_to_graph Length FRC   SpeedLimit StreetName     
##    <int> <dbl> <dbl>              <dbl>  <dbl> <chr>      <dbl> <chr>          
##  1    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  2    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  3    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  4  6666   0       1            0.00100 0.0361 5             40 Harrison St    
##  5    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  6    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  7    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  8    NA  NA      NA           NA       0.0361 5             40 Harrison St    
##  9    NA  NA      NA           NA       0.0361 5             40 Harrison St    
## 10  8941  12.9     1            0.00230 0.112  6             35 Rhode Island St
## # ℹ 57,098 more rows
## # ℹ 46 more variables: harmonicAverageSpeed <dbl>, medianSpeed <dbl>,
## #   averageSpeed <dbl>, sampleSize <int>, averageTravelTime <dbl>,
## #   medianTravelTime <dbl>, travelTimeRatio <dbl>, List_Number <int>,
## #   `5percentile` <int>, `10percentile` <int>, `15percentile` <int>,
## #   `20percentile` <int>, `25percentile` <int>, `30percentile` <int>,
## #   `35percentile` <int>, `40percentile` <int>, `45percentile` <int>, …
summary(sf_graph)
## A metric graph object with:
## 
## Vertices:
##   Total: 4017 
##   Degree 2: 1821;  Degree 3: 180;  Degree 4: 1402;  Degree 5: 94;  Degree 6: 367; 
##   Degree 7: 36;  Degree 8: 116;  Degree 12: 1; 
##   With incompatible directions:  0 
## 
## Edges: 
##   Total: 6827 
##   Lengths: 
##       Min: 0.002834658  ; Max: 0.2743201  ; Total: 311.1441 
##   Weights: 
##       Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour 
##   That are circles:  0 
## 
## Graph units: 
##   Vertices unit:  degrees  ; Lengths unit:  km 
## 
## Longitude and Latitude coordinates:  TRUE
##   Which spatial package:  sf 
##   CRS:  EPSG:4326
## 
## Some characteristics of the graph:
##   Connected: TRUE
##   Has loops: FALSE
##   Has multiple edges: TRUE
##   Is a tree: FALSE
##   Distance consistent: TRUE
##   Has Euclidean edges: FALSE
## 
## Computed quantities inside the graph: 
##   Laplacian:  FALSE  ; Geodesic distances:  TRUE 
##   Resistance distances:  FALSE  ; Finite element matrices:  FALSE 
## 
## Mesh: The graph has no mesh! 
## 
## Data: 
##   Columns:  ID speed day Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour 
##   Groups:  .group 
## 
## Tolerances: 
##   vertex-vertex:  0.001 
##   vertex-edge:  0.001 
##   edge-edge:  0

When running sf_graph$edgeweight_to_data(data_loc = TRUE), some NA values are created (because the data is grouped). We remove them below. We also standardize the SpeedLimit covariate.

data = sf_graph$get_data() %>% 
  drop_na(-StreetName) %>% # this drops all rows with at least one NA value but without taking into account StreetName
  mutate(across(c("SpeedLimit"), ~standardize(.))) %>%
  dplyr::select(speed, SpeedLimit)

The code of chunk below was executed only one time.


{r, eval = FALSE}
aux = data |>
  rename(distance_on_edge = .distance_on_edge, edge_number = .edge_number) |>
  as.data.frame() |>
  dplyr::select(edge_number, distance_on_edge, .group)

distmatrixlist = list()

for (i in 1:4) {
  distmatrixlist[[i]] = sf_graph$compute_geodist_PtE(PtE = aux %>% 
                                                       filter(.group == as.character(i)) %>% 
                                                       dplyr::select(-.group),
                                                     normalized = TRUE,
                                                     include_vertices = FALSE)
}


save(distmatrixlist, file = here("Models_output/distmatrixfixed.RData"))

The code of chunk above was executed only one time.


summary(data)
##      speed          SpeedLimit         .group           .edge_number 
##  Min.   : 0.000   Min.   :-2.3743   Length:14534       Min.   :   1  
##  1st Qu.: 9.656   1st Qu.:-0.1006   Class :character   1st Qu.:1376  
##  Median :19.312   Median :-0.1006   Mode  :character   Median :2956  
##  Mean   :19.385   Mean   : 0.0000                      Mean   :3130  
##  3rd Qu.:27.359   3rd Qu.:-0.1006                      3rd Qu.:4803  
##  Max.   :99.779   Max.   : 6.6173                      Max.   :6817  
##  .distance_on_edge    .coord_x         .coord_y    
##  Min.   :0.0000    Min.   :-122.4   Min.   :37.77  
##  1st Qu.:0.2552    1st Qu.:-122.4   1st Qu.:37.78  
##  Median :0.5139    Median :-122.4   Median :37.79  
##  Mean   :0.5099    Mean   :-122.4   Mean   :37.79  
##  3rd Qu.:0.7697    3rd Qu.:-122.4   3rd Qu.:37.79  
##  Max.   :1.0000    Max.   :-122.4   Max.   :37.81

We add the data again but now with the new standardized SpeedLimit covariate.

sf_graph$add_observations(data = data, 
                          group = "day", 
                          normalized = TRUE, 
                          clear_obs = TRUE)
sf_graph$get_data()
## # A tibble: 14,534 × 7
##    speed SpeedLimit .coord_x .coord_y .edge_number .distance_on_edge .group
##    <dbl>      <dbl>    <dbl>    <dbl>        <dbl>             <dbl> <chr> 
##  1   0       -0.101    -122.     37.8            1            0.437  1     
##  2  12.9     -0.617    -122.     37.8            4            0.144  1     
##  3  24.1     -0.617    -122.     37.8            6            0.252  1     
##  4  32.2     -0.617    -122.     37.8            6            0.658  1     
##  5   0       -0.927    -122.     37.8            9            0.601  1     
##  6  19.3     -0.927    -122.     37.8           14            0.0247 1     
##  7  22.5     -0.927    -122.     37.8           14            0.362  1     
##  8  19.3     -0.927    -122.     37.8           14            0.832  1     
##  9  30.6     -0.927    -122.     37.8           18            0.358  1     
## 10  14.5     -0.927    -122.     37.8           20            0.309  1     
## # ℹ 14,524 more rows
summary(sf_graph)
## A metric graph object with:
## 
## Vertices:
##   Total: 4017 
##   Degree 2: 1821;  Degree 3: 180;  Degree 4: 1402;  Degree 5: 94;  Degree 6: 367; 
##   Degree 7: 36;  Degree 8: 116;  Degree 12: 1; 
##   With incompatible directions:  0 
## 
## Edges: 
##   Total: 6827 
##   Lengths: 
##       Min: 0.002834658  ; Max: 0.2743201  ; Total: 311.1441 
##   Weights: 
##       Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour 
##   That are circles:  0 
## 
## Graph units: 
##   Vertices unit:  degrees  ; Lengths unit:  km 
## 
## Longitude and Latitude coordinates:  TRUE
##   Which spatial package:  sf 
##   CRS:  EPSG:4326
## 
## Some characteristics of the graph:
##   Connected: TRUE
##   Has loops: FALSE
##   Has multiple edges: TRUE
##   Is a tree: FALSE
##   Distance consistent: TRUE
##   Has Euclidean edges: FALSE
## 
## Computed quantities inside the graph: 
##   Laplacian:  FALSE  ; Geodesic distances:  TRUE 
##   Resistance distances:  FALSE  ; Finite element matrices:  FALSE 
## 
## Mesh: The graph has no mesh! 
## 
## Data: 
##   Columns:  speed SpeedLimit 
##   Groups:  .group 
## 
## Tolerances: 
##   vertex-vertex:  0.001 
##   vertex-edge:  0.001 
##   edge-edge:  0

We build a mesh.

h = 0.05
sf_graph$build_mesh(h = h)
summary(sf_graph)
## A metric graph object with:
## 
## Vertices:
##   Total: 4017 
##   Degree 2: 1821;  Degree 3: 180;  Degree 4: 1402;  Degree 5: 94;  Degree 6: 367; 
##   Degree 7: 36;  Degree 8: 116;  Degree 12: 1; 
##   With incompatible directions:  0 
## 
## Edges: 
##   Total: 6827 
##   Lengths: 
##       Min: 0.002834658  ; Max: 0.2743201  ; Total: 311.1441 
##   Weights: 
##       Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour 
##   That are circles:  0 
## 
## Graph units: 
##   Vertices unit:  degrees  ; Lengths unit:  km 
## 
## Longitude and Latitude coordinates:  TRUE
##   Which spatial package:  sf 
##   CRS:  EPSG:4326
## 
## Some characteristics of the graph:
##   Connected: TRUE
##   Has loops: FALSE
##   Has multiple edges: TRUE
##   Is a tree: FALSE
##   Distance consistent: TRUE
##   Has Euclidean edges: FALSE
## 
## Computed quantities inside the graph: 
##   Laplacian:  FALSE  ; Geodesic distances:  TRUE 
##   Resistance distances:  FALSE  ; Finite element matrices:  FALSE 
## 
## Mesh: 
##   Max h_e:  0.04999798  ; Min n_e:  0 
## 
## Data: 
##   Columns:  speed SpeedLimit 
##   Groups:  .group 
## 
## Tolerances: 
##   vertex-vertex:  0.001 
##   vertex-edge:  0.001 
##   edge-edge:  0

We get the value of the weights at mesh locations. This will allow us to built matrices B.sigma and B.range below. Again, sf_graph$edgeweight_to_data(mesh = TRUE, add = FALSE, return = TRUE) creates repeated information (because the data is grouped). We fix that by filtering one group. We also standardize the SpeedLimit covariate.

mesh = sf_graph$edgeweight_to_data(mesh = TRUE, 
                                   add = FALSE, 
                                   return = TRUE) %>% 
  filter(.group == 1) %>%
  mutate(across(c("SpeedLimit"), ~standardize(.))) %>%
  dplyr:::select.data.frame(SpeedLimit)
summary(mesh)
##    SpeedLimit     
##  Min.   :-1.9800  
##  1st Qu.:-0.1744  
##  Median :-0.1744  
##  Mean   : 0.0000  
##  3rd Qu.:-0.1744  
##  Max.   : 5.1604

1.1 Stationary model

  • Observe that we are considering replicates.
stat.time.ini <- Sys.time()
################################################################################
################################# STATIONARY MODEL #############################
################################################################################

rspde_model_stat <- rspde.metric_graph(sf_graph,
                                         parameterization = "matern",
                                         nu.upper.bound = 1.5)
str(rspde_model_stat)
## List of 20
##  $ f                   :List of 3
##   ..$ model   : chr "cgeneric"
##   ..$ n       : int 21507
##   ..$ cgeneric:List of 5
##   .. ..$ model: chr "inla_cgeneric_rspde_stat_general_model"
##   .. ..$ shlib: chr "/usr/local/lib/R/site-library/INLA/bin/linux/64bit/external/rSPDE/librSPDE.so"
##   .. ..$ n    : int 21507
##   .. ..$ debug: logi FALSE
##   .. ..$ data :List of 5
##   .. .. ..$ ints      :List of 5
##   .. .. .. ..$ n          : int 21507
##   .. .. .. ..$ debug      : int 0
##   .. .. .. ..$ graph_opt_i: int [1:135512] 0 0 0 0 0 0 0 0 0 0 ...
##   .. .. .. ..$ graph_opt_j: int [1:135512] 0 1 230 984 985 986 1700 1701 2280 2455 ...
##   .. .. .. ..$ rspde.order: int 2
##   .. .. ..$ doubles   :List of 11
##   .. .. .. ..$ d                   : num 1
##   .. .. .. ..$ nu.upper.bound      : num 1.5
##   .. .. .. ..$ matrices_less       : num [1:94764] 0.0797 0 0 0 0 ...
##   .. .. .. ..$ matrices_full       : num [1:207848] 0.0797 0 0 0 0 ...
##   .. .. .. ..$ start.theta         : num [1:2] 0 0.223
##   .. .. .. ..$ theta.prior.mean    : num [1:2] 0 0.223
##   .. .. .. ..$ prior.nu.loglocation: num -0.288
##   .. .. .. ..$ prior.nu.mean       : num 0.75
##   .. .. .. ..$ prior.nu.prec       : num 3
##   .. .. .. ..$ prior.nu.logscale   : num 1
##   .. .. .. ..$ start.nu            : num 0.75
##   .. .. ..$ characters:List of 5
##   .. .. .. ..$ model            : chr "inla_cgeneric_rspde_stat_general_model"
##   .. .. .. ..$ shlib            : chr "/usr/local/lib/R/site-library/INLA/bin/linux/64bit/external/rSPDE/librSPDE.so"
##   .. .. .. ..$ prior.nu.dist    : chr "lognormal"
##   .. .. .. ..$ parameterization : chr "matern"
##   .. .. .. ..$ prior.theta.param: chr "theta"
##   .. .. ..$ matrices  :List of 2
##   .. .. .. ..$ rational_table  : num [1:5996] 999 6 0.001 0.412 0.005 ...
##   .. .. .. ..$ theta.prior.prec: num [1:6] 2 2 0.1 0 0 0.1
##   .. .. ..$ smatrices : list()
##   .. ..- attr(*, "class")= chr "inla.cgeneric"
##  $ cgeneric_type       : chr "general"
##  $ theta.prior.mean    : num [1:2] 0 0.223
##  $ prior.nu            :List of 4
##   ..$ loglocation: num -0.288
##   ..$ mean       : num 0.75
##   ..$ prec       : num 3
##   ..$ logscale   : num 1
##  $ theta.prior.prec    : num [1:2, 1:2] 0.1 0 0 0.1
##  $ start.nu            : num 0.75
##  $ integer.nu          : logi FALSE
##  $ start.theta         : num [1:2] 0 0.223
##  $ stationary          : logi TRUE
##  $ rspde.order         : num 2
##  $ dim                 : num 1
##  $ est_nu              : logi TRUE
##  $ nu.upper.bound      : num 1.5
##  $ prior.nu.dist       : chr "lognormal"
##  $ debug               : logi FALSE
##  $ type.rational.approx: chr "chebfun"
##  $ mesh                :Classes 'metric_graph', 'R6' <metric_graph>
##   Public:
##     add_mesh_observations: function (data = NULL, group = NULL) 
##     add_observations: function (data = NULL, edge_number = "edge_number", distance_on_edge = "distance_on_edge", 
##     build_mesh: function (h = NULL, n = NULL, continuous = TRUE, continuous.outs = FALSE, 
##     buildC: function (alpha = 2, edge_constraint = FALSE) 
##     buildDirectionalConstraints: function (alpha = 1) 
##     C: NULL
##     characteristics: list
##     check_distance_consistency: function () 
##     check_euclidean: function () 
##     clear_observations: function () 
##     clone: function (deep = FALSE) 
##     CoB: NULL
##     compute_characteristics: function (check_euclidean = FALSE) 
##     compute_fem: function (petrov = FALSE) 
##     compute_geodist: function (full = FALSE, obs = TRUE, group = NULL, verbose = 0) 
##     compute_geodist_mesh: function () 
##     compute_geodist_PtE: function (PtE, normalized = TRUE, include_vertices = TRUE, verbose = 0) 
##     compute_laplacian: function (full = FALSE, obs = TRUE, group = NULL, verbose = 0) 
##     compute_PtE_edges: function () 
##     compute_resdist: function (full = FALSE, obs = TRUE, group = NULL, check_euclidean = FALSE, 
##     compute_resdist_mesh: function () 
##     compute_resdist_PtE: function (PtE, normalized = TRUE, include_vertices = FALSE, check_euclidean = FALSE, 
##     coordinates: function (PtE = NULL, XY = NULL, normalized = TRUE) 
##     drop_na: function (...) 
##     E: 1 3 3 6 4 7 9 11 13 14 15 16 18 19 20 21 23 25 27 29 31  ...
##     edge_lengths: 0.0363234139144278 0.01586106867077 0.027923568765887 0. ...
##     edges: metric_graph_edges
##     edgeweight_to_data: function (loc = NULL, mesh = FALSE, data_loc = FALSE, weight_col = NULL, 
##     fem_basis: function (PtE) 
##     filter: function (..., .drop_na = FALSE, .drop_all_na = TRUE) 
##     geo_dist: list
##     get_data: function (group = NULL, tibble = TRUE, drop_na = FALSE, drop_all_na = TRUE) 
##     get_degrees: function (which = "degree") 
##     get_edge_lengths: function (unit = NULL) 
##     get_edge_weights: function (data.frame = FALSE, tibble = TRUE) 
##     get_groups: function (get_cols = FALSE) 
##     get_initial_graph: function () 
##     get_locations: function () 
##     get_mesh_locations: function (bru = FALSE, loc = NULL, normalized = TRUE) 
##     get_PtE: function () 
##     get_vertices_incomp_dir: function () 
##     initialize: function (edges = NULL, V = NULL, E = NULL, vertex_unit = NULL, 
##     is_tree: function () 
##     Laplacian: NULL
##     mesh: list
##     mesh_A: function (PtE) 
##     mutate: function (..., .drop_na = FALSE, .drop_all_na = TRUE) 
##     nE: 6827
##     nV: 4017
##     observation_to_vertex: function (tolerance = 1e-15, mesh_warning = TRUE) 
##     plot: function (data = NULL, newdata = NULL, group = 1, plotly = FALSE, 
##     plot_connections: function () 
##     plot_function: function (data = NULL, newdata = NULL, group = 1, X = NULL, plotly = FALSE, 
##     plot_movie: function (X, plotly = TRUE, vertex_size = 5, vertex_color = "black", 
##     print: function () 
##     process_data: function (data = NULL, edge_number = "edge_number", distance_on_edge = "distance_on_edge", 
##     prune_vertices: function (check_weights = TRUE, verbose = FALSE) 
##     PtV: NULL
##     res_dist: NULL
##     select: function (..., .drop_na = FALSE, .drop_all_na = TRUE) 
##     set_edge_weights: function (weights = rep(1, self$nE), kirchhoff_weights = NULL) 
##     summarise: function (..., .include_graph_groups = FALSE, .groups = NULL, 
##     summary: function (messages = FALSE, compute_characteristics = TRUE, check_euclidean = TRUE, 
##     V: -122.41277 -122.41249 -122.40376 -122.40358 -122.40379 - ...
##     vertices: metric_graph_vertices
##     VtEfirst: function () 
##   Private:
##     A: function (group = NULL, obs_to_vert = FALSE, drop_na = FALSE, 
##     add_vertices: function (PtE, tolerance = 1e-10, verbose) 
##     addinfo: FALSE
##     clear_initial_info: function () 
##     compute_degrees: function () 
##     compute_laplacian_PtE: function (PtE, normalized = TRUE, verbose = verbose) 
##     compute_lengths: function (longlat, unit, crs, proj4string, which_longlat, vertex_unit, 
##     connected: TRUE
##     coordinates_multiple_snaps: function (XY, tolerance, verbose = verbose, crs, proj4string, 
##     create_update_vertices: function () 
##     crs: crs
##     data: metric_graph_data, list
##     edge_weights: tbl_df, tbl, data.frame
##     find_edge_edge_points: function (tol, verbose, crs, proj4string, longlat, fact, which_longlat) 
##     find_mesh_bc: function () 
##     get_edge_weights_internal: function (data.frame = FALSE) 
##     group_col: .group
##     initial_edges_added: NULL
##     initial_graph: metric_graph, R6
##     kirchhoff_weights: NULL
##     length_unit: km
##     line_to_vertex: function (tolerance = 0, longlat = FALSE, fact, verbose, crs, 
##     longlat: TRUE
##     merge_close_vertices: function (tolerance, fact) 
##     merge.all.deg2: function () 
##     mesh_merge_deg2: function () 
##     mesh_merge_outs: function () 
##     move_V_first: function () 
##     plot_2d: function (line_width = 0.1, marker_size = 1, vertex_color = "black", 
##     plot_3d: function (line_width = 1, marker_size = 1, vertex_color = "rgb(0,0,0)", 
##     proj4string: NULL
##     prune_warning: FALSE
##     pruned: FALSE
##     PtE_to_mesh: function (PtE) 
##     ref_edges: 1 1552 2 5 3133 4 6 395 7 526 8 6466 9 10 11 12 38 13 14 ...
##     remove_circles: function (threshold, verbose, longlat, unit, crs, proj4string, 
##     remove.first.deg2: function (res) 
##     set_first_weights: function (weights = rep(1, self$nE)) 
##     set_petrov_matrices: function () 
##     split_edge: function (Ei, t, tolerance = 0) 
##     temp_PtE: NULL
##     tolerance: list
##     transform: FALSE
##     vertex_unit: degrees
##     which_longlat: sf 
##  $ fem_mesh            :List of 5
##   ..$ c0:Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
##   .. .. ..@ i       : int [1:7169] 0 1 2 3 4 5 6 7 8 9 ...
##   .. .. ..@ j       : int [1:7169] 0 1 2 3 4 5 6 7 8 9 ...
##   .. .. ..@ Dim     : int [1:2] 7169 7169
##   .. .. ..@ Dimnames:List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : NULL
##   .. .. ..@ x       : num [1:7169] 0.0797 0.1564 0.1117 0.0493 0.0614 ...
##   .. .. ..@ factors : list()
##   ..$ g1:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. ..@ i       : int [1:24879] 0 1 985 0 1 1700 5207 5364 6858 7041 ...
##   .. .. ..@ p       : int [1:7170] 0 3 10 16 20 23 27 36 44 50 ...
##   .. .. ..@ Dim     : int [1:2] 7169 7169
##   .. .. ..@ Dimnames:List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : NULL
##   .. .. ..@ x       : num [1:24879] 101.1 -55.1 -46.1 -55.1 216.4 ...
##   .. .. ..@ factors : list()
##   ..$ g2:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. ..@ i       : int [1:56007] 0 1 985 986 1700 4722 5207 5364 6240 6858 ...
##   .. .. ..@ p       : int [1:7170] 0 11 22 34 43 52 61 76 89 99 ...
##   .. .. ..@ Dim     : int [1:2] 7169 7169
##   .. .. ..@ Dimnames:List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : NULL
##   .. .. ..@ x       : num [1:56007] 172097 -145986 -232979 119006 15770 ...
##   .. .. ..@ factors : list()
##   ..$ g3:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. ..@ i       : int [1:96755] 0 1 230 984 985 986 1700 1701 2280 2455 ...
##   .. .. ..@ p       : int [1:7170] 0 17 35 57 76 91 103 126 150 166 ...
##   .. .. ..@ Dim     : int [1:2] 7169 7169
##   .. .. ..@ Dimnames:List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : NULL
##   .. .. ..@ x       : num [1:96755] 3.94e+08 -3.74e+08 -2.65e+07 -9.14e+07 -1.47e+09 ...
##   .. .. ..@ factors : list()
##   ..$ g4:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. ..@ i       : int [1:147901] 0 1 230 983 984 985 986 989 1700 1701 ...
##   .. .. ..@ p       : int [1:7170] 0 28 52 89 122 151 171 209 249 274 ...
##   .. .. ..@ Dim     : int [1:2] 7169 7169
##   .. .. ..@ Dimnames:List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : NULL
##   .. .. ..@ x       : num [1:147901] 1.41e+12 -1.06e+12 -2.25e+11 1.91e+11 -1.41e+12 ...
##   .. .. ..@ factors : list()
##  $ parameterization    : chr "matern"
##  $ n.spde              : int 7169
##  - attr(*, "class")= chr [1:3] "rspde_metric_graph" "inla_rspde" "inla.cgeneric"
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
                                        repl = ".all",
                                        loc_name = "loc")
str(data_rspde_bru_stat)
## List of 4
##  $ data :List of 8
##   ..$ speed            : num [1:14534] 0 12.9 24.1 32.2 0 ...
##   ..$ SpeedLimit       : num [1:14534] -0.101 -0.617 -0.617 -0.617 -0.927 ...
##   ..$ .coord_x         : num [1:14534] -122 -122 -122 -122 -122 ...
##   ..$ .coord_y         : num [1:14534] 37.8 37.8 37.8 37.8 37.8 ...
##   ..$ .edge_number     : num [1:14534] 1 4 6 6 9 14 14 14 18 20 ...
##   ..$ .distance_on_edge: num [1:14534] 0.437 0.144 0.252 0.658 0.601 ...
##   ..$ .group           : chr [1:14534] "1" "1" "1" "1" ...
##   ..$ loc              : num [1:14534, 1:2] 1 4 6 6 9 14 14 14 18 20 ...
##   ..- attr(*, "class")= chr [1:2] "metric_graph_data" "list"
##  $ index:List of 3
##   ..$ field      : int [1:28676] 1 2 3 4 5 6 7 8 9 10 ...
##   ..$ field.group: int [1:28676] 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ field.repl : int [1:28676] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "class")= chr [1:2] "inla_rspde_index" "list"
##   ..- attr(*, "rspde.order")= num 0
##   ..- attr(*, "integer_nu")= logi TRUE
##   ..- attr(*, "n.mesh")= int 7169
##   ..- attr(*, "name")= chr "field"
##   ..- attr(*, "n.group")= int 1
##   ..- attr(*, "n.repl")= int 4
##  $ repl : chr [1:14534] "1" "1" "1" "1" ...
##  $ basis:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. ..@ i       : int [1:87204] 0 555 1905 1906 0 1050 1471 1472 1473 1905 ...
##   .. ..@ p       : int [1:86029] 0 4 11 11 11 11 12 13 15 15 ...
##   .. ..@ Dim     : int [1:2] 14534 86028
##   .. ..@ Dimnames:List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : NULL
##   .. ..@ x       : num [1:87204] 0.563 0.0495 0.3595 0.7935 0.437 ...
##   .. ..@ factors : list()
cmp_stat = speed ~ -1 +
  Intercept(1) +
  SpeedLimit +
  field(loc, model = rspde_model_stat,
        replicate = data_rspde_bru_stat[["repl"]])

rspde_fit_stat <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )
## Error in inla.inlaprogram.has.crashed() : 
##   The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
##   If this does not help, please contact the developers at <help@r-inla.org>.
## 
##  *** inla.core.safe:  inla.program has crashed: rerun to get better initial values. try=1/2 
## Error in inla.inlaprogram.has.crashed() : 
##   The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
##   If this does not help, please contact the developers at <help@r-inla.org>.
## 
##  *** inla.core.safe:  inla.program has crashed: rerun to get better initial values. try=2/2 
## Error in inla.inlaprogram.has.crashed() : 
##   The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
##   If this does not help, please contact the developers at <help@r-inla.org>.
str(rspde_fit_stat)
## Error in eval(expr, envir, enclos): object 'rspde_fit_stat' not found
stat.time.fin <- Sys.time()
print(stat.time.fin - stat.time.ini)
## Time difference of 1.323035 mins
summary(rspde_fit_stat)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rspde_fit_stat' not found
fit.rspde = rspde.result(rspde_fit_stat, "field", rspde_model_stat)
## Error in eval(expr, envir, enclos): object 'rspde_fit_stat' not found
summary(fit.rspde)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fit.rspde' not found

1.2 Nonstationary model

  • Observe that we are using the computed parameters from the stationary model as initial values for the nonstationary models.
nonstat.time.ini <- Sys.time()
################################################################################
############################# NON STATIONARY MODEL #############################
################################################################################

B.sigma = cbind(0, 1, 0, mesh$SpeedLimit, 0)
B.range = cbind(0, 0, 1, 0, mesh$SpeedLimit)
init.vec.theta = c(fit.rspde$summary.log.std.dev$mode, 
                   fit.rspde$summary.log.range$mode, 
                   rep(0, (ncol(B.sigma)-3)))
## Error in eval(expr, envir, enclos): object 'fit.rspde' not found

rspde_model_nonstat <- rspde.metric_graph(sf_graph,
                                          start.theta = init.vec.theta,
                                          theta.prior.mean = init.vec.theta,
                                          B.sigma = B.sigma,
                                          B.range = B.range,
                                          parameterization = "matern",
                                          nu.upper.bound = 1.5)
## Error in eval(expr, envir, enclos): object 'init.vec.theta' not found
str(rspde_model_nonstat)
## Error in eval(expr, envir, enclos): object 'rspde_model_nonstat' not found
data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
                                           repl = ".all",
                                           loc_name = "loc")
## Error in eval(expr, envir, enclos): object 'rspde_model_nonstat' not found
str(data_rspde_bru_nonstat)
## Error in eval(expr, envir, enclos): object 'data_rspde_bru_nonstat' not found
cmp_nonstat = speed ~ -1 +
  Intercept(1) +
  SpeedLimit +
  field(loc, model = rspde_model_nonstat,
        replicate = data_rspde_bru_nonstat[["repl"]])

rspde_fit_nonstat <-
  bru(cmp_nonstat,
      data = data_rspde_bru_nonstat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )
## Error in eval(expr, envir, enclos): object 'data_rspde_bru_nonstat' not found
str(rspde_fit_nonstat)
## Error in eval(expr, envir, enclos): object 'rspde_fit_nonstat' not found
nonstat.time.fin <- Sys.time()
print(nonstat.time.fin - nonstat.time.ini)
## Time difference of 0.02544212 secs
summary(rspde_fit_nonstat)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rspde_fit_nonstat' not found
summary(rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rspde_model_nonstat' not found

1.3 Crossvalidation 1

#load(here("Models_output/distmatrixfixed.RData"))

points = data %>%
  as.data.frame() %>%
  st_as_sf(coords = c(".coord_x", ".coord_y"), crs = 4326) %>%
  mutate(., index = 1:nrow(.)) %>% 
  st_drop_geometry() %>%
  dplyr:::select(speed, .group, index) %>%
  mutate(.group = as.numeric(.group)) %>%
  group_by(.group) %>%
  mutate(indexingroup = seq_len(n())) %>%
  ungroup()

distance = seq(from = 0, to = 400, by = 20)/1000

#save(GROUPS, file = here("Models_output/GROUPS_for_window_case.RData"))
load(here("Models_output/GROUPS_for_window_case.RData"))
mse.stat <- mse.nonstat <- ls.stat <- ls.nonstat <- rep(0,length(distance))
# cross-validation for-loop
for (j in 1:length(distance)) {
  print(j)
  # cross-validation of the stationary model
  cv.stat <- inla.group.cv(rspde_fit_stat, groups = GROUPS[[j]])
  # cross-validation of the nonstationary model
  cv.nonstat <- inla.group.cv(rspde_fit_nonstat, groups = GROUPS[[j]])
  # obtain MSE and LS
  mse.stat[j] <- mean((cv.stat$mean - points$speed)^2)
  mse.nonstat[j] <- mean((cv.nonstat$mean - points$speed)^2)
  ls.stat[j] <- mean(log(cv.stat$cv))
  ls.nonstat[j] <- mean(log(cv.nonstat$cv))
}
## [1] 1
## Error in eval(expr, envir, enclos): object 'rspde_fit_stat' not found

## plot results
par(mfrow = c(2,2), family = "Palatino")

# Plot MSE
plot(distance, mse.stat, main = "MSE", ylim = c(min(mse.nonstat, mse.stat), max(mse.nonstat, mse.stat)),
     type = "l", ylab = "MSE", xlab = "distance in m", col = "black")
lines(distance, mse.nonstat, col = "blue")
legend("bottomright", legend = c("Stationary", "Non-stationary"), col = c("black", "blue"), lty = 1)

# Plot log-score
plot(distance, -ls.stat, main = "log-score", ylim = c(min(-ls.nonstat, -ls.stat), max(-ls.nonstat, -ls.stat)),
     type = "l", ylab = "log-score", xlab = "distance in m", col = "black")
lines(distance, -ls.nonstat, col = "blue")
legend("bottomright", legend = c("Stationary", "Non-stationary"), col = c("black", "blue"), lty = 1)

save.image(here(paste0("Models_output/", rmarkdown::metadata$title, ".RData")))